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§1.  STATEMENT  OF  SCIENTIFIC  WORK  DURING  THE  REPORTING  °ER10D. 

Let  ”  be  the  N  dimensional  real  euclidean  space, 
/  T  r>xl 

x^  =  (x1  ,  x^  ,  ..  ,  X^)  fe K  a  generic  vector,  the  superscript 

T  means  transpose,  (■  ,  )  is  the  euclidean  scalar  product  and 

(\  *  l\  the  euclidean  norm. 

We  have  considered  the  linear  programming  problem; 
that  i s : 

Problem  1  (Linear  programming)  Given  b,  c  fc  IR  and  A  an  mxhl 
matrix  ( m  <  N)  ,  solve  the  following  minimization  problem: 


minimize  <  c  ,  _»> 
for  x  such  that 

X  i  0 


Ax  -  b  ^  0 


where  x^  >  Q  means  that  each  component  of  x  is  greater  or 
equal  to  zero  and  similarly  Ax  -  b  0. 

The  Linear  programming  problem  is  transformed  into  a 
Linear  Complementarity  Problem  using  the  Lagrange  multipliers 
as  follows: 


i  -£  \ 


-aT  x 


Problem  1  becomes 
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Problem  2  (Linear  Complementarity  Problem).  Find  z^  such  that: 


L  Z-V 

/N.  A/  . 

A2  -  t  j  o 
(  i*  Az  -  b  >  =  0 


Finally  the  Linear  Complementarity  Problem  is  trans¬ 
formed  into  a  global  optimization  problem  via  the  transforma¬ 
tion  given  at  page  1  of  the  Fourth  Periodic  Report. In  this  way  the  li¬ 
near  programming  problem  is  reduced  to  the  problem  of  finding 
the  minimizers  of  a  function  F  (  z  )  '•?•  0  such  that  F  =  0  at  the 
minimizers  . 

If  the  original  linear  programming  problem  has  a 
unique  solution  the  corresponding  minimization  problem  has 
a  unique  global  minimizer,  but  unfortunately  the  function  F 
may  have  a  rather  complicated  set  of  local  minimizers  inclu¬ 
ding  "cilindrical  valleys".  A  direct  application  of  a  local 
minimization  techniques  such  as  GRACON  (see  Fourth  Periodic 
Report)  is  unsuccessful.  So  that  several  modified  algorithms 
have  been  considered  the  most  promising  ones  are: 

(i)  A  perturbation  method.  Let  £ e  [R.  ^ £>£  ,  the  matrix  A  is 
substituted  with  A  (£.)  =  A  +  £l  where  1  is  the  identity  ma¬ 
trix.  The  idea  is  to  find  the  solution  zit)  of  the  prob- 
lem  correspondent  to  the  matrix  A(£)  and  to  let  £.  go  to 
zero.  If  £  goes  to  zero  too  quickly  the  conjugate  gra¬ 
dient  method  may  remain  trapped  in  a  local  minimizer,  if 
£  goes  to  zero  too  slowly  the  computational  cost  becomes 
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excessive.  In  Table  1  the  results  obtained  with  this  me¬ 
thod  are  indicated  with  00D4  DOD; 

(ii)  Avoiding  the  cilindrical  valleys.  When  the  conjugate  gra 
dient  method  remains  trapped  in  a  cilindrical  valley  the 
minimization  procedure  is  stopped,  the  direction  of  the 
valley  is  computed  and  a  large  step  is  taken  in  the  di¬ 
rection  of  the  valley  in  the  "descending"  direction,  fi¬ 
nally  the  conjugate  gradient  minimization  algorithm  is 
started  again.  In  Table  1  the  results  obtained  with  this 
method  are  indicated  with  00D4  D009. 

The  test  problems  considered  in  Table  1  have  been  con¬ 
structed  following  the  suggestions  of  De  Leone  and  Manga- 
sarian  (see  Appendix  2).  The  numerical  experience  has  been 
obtained  on  a  VAX  8530  with  VMS  4,5  Operating  System. 
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TABLE  1 


0004 

0CD 

0004 

D009 

M 

NFEV 

T 

M 

NFEV 

T 

50 

694 

2" 

.  97 

50 

315 

2" 

.  24 

100 

1502 

7" 

.87 

100 

716 

4" 

.21 

500 

4176 

1  '  .35" 

.  08 

500 

2154 

43" 

.  44 

1000 

9620 

7'.  4" 

.74 

1000 

1333 

55" 

.  39 

5000 

1  2548 

46  '  .  26" 

.80 

10000 

14241 

1h. 54 ' . 16" 

.  51 

Legenda 

M  dimension  of  the  complementarity  problem  (z  £  IK  ) 

NFEV  number  of  function  and  gradient  evaluation  of  the  func¬ 
tion  that  must  be  minimized 

T  running  time  on  the  VAX  8530  (Example  I*1  .54'. 18". 51  =  1 

hour  54  minutes,  18  seconds,  51  seconds/100. 

Tolerance  required:  Final  function  value  ^  10 


Finally  some  geophysical  application  of  linear  pro¬ 
gramming  has  been  considered  (see  Appendix  1)  since  in  our 
opinion  the  random  generated  test  problems  are  not  completely 
satisfactory  . 


§2.  RESEARCH  PLANS  FOR  THE  IMMEDIATE  FUTURE. 


In  the  immediate  future  we  plan  to  pursue  the  fol 
lowing  objectives: 

(i)  study  Karmarkar  and  Renegar  methods  for  linear  program 
ming  in  the  context  of  continuation  methods; 

( i i )  study  the  nonlinear  complementarity  problem. 
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§3.  ADMINISTRATIVE  ACTIONS. 

None  . 
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§4.  Appendix  1:  L.  Misici,  F.  Zirilli:  "The  inverse  gravi¬ 
metry  problem:  An  application  to  the  northern  San  Fran¬ 
cisco  craton  granite"  submitted  to  J.  of  the  Geological 
Society  London. 
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THE  INVERSE  GRAVIMETRY  PROBLEM:  AN  APPLICATION 
TO  THE  NORTHERN  SAN  FRANCISCO  CRATON  GRANITE* 


Luciano  Misici 

Dipartimento  ii  Matematica  e  Fistca 
Universita  ii  Camerino 
SSOSS  Camerino  (MC)  Italy 


Francesco  ZLrilli 

Dipartimento  ii  Matematica  ’G.  Castelnuovo ’ 
Universita  ii  Roma  ’La  Sapienza’ 

00185  Roma  Italy 


Abstract 

From  the  knowledge  of  the  anomalies  of  the  gravitational  field  we  reconstruct  the  mass  density 
distribution  in  a  large  region  of  the  state  of  Bahia  (Brazil). This  inverse  gravimetry  problem  has 
been  translated  in  a  linear  programming  problem  and  solved  using  the  simplex  method.  Both  two 
and  three  dimensional  models  have  been  considered. 

1  •  Introduction 

In  recent  years  a  great  deal  of  attention  has  been  attracted  by  the  study  of  several  inverse 
problems  in  science  and  technology.  In  particular  in  science  we  mention  the  inverse  problem  of 
quantum  mechanics,  that  is  the  reconstruction  of  a  potential  from  its  scattering  data  [lj  and  the 
famous  problem  "can  you  hear  the  shape  of  a  drum”,  that  is  the  reconstruction  of  the  shape 
of  a  region  in  Euclidean  space  from  the  knowledge  of  spectral  properties  of  suitable  differential 
operators  (2j.  In  technology  we  mention  the  radar  technology  that  is  the  use  of  electromagnetic 
waves  to  detect  not  visible  objects,  the  sonar  and  ecography  technology  that  is  the  use  of  acoustics 
waves  to  detect  properties  of  regions  that  are  not  directly  accessible  because  they  are  underwater 
or  inside  the  human  body. 

*  The  research  reported  in  this  document  has  been  made  possible  through  the  support  and  sponsorship  of  the 
U.S.  Government  through  the  European  Research  Office  of  the  U.S.  Array  under  contract  n.  DAJA  45-8C-C-002S 
at  the  University  of  Roma  "La  Sapienta"  and  of  the  Ministero  Pubblica  Istruiione  under  contract  G0°%  1987  at  the 
University  of  Camerino. 
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There  ire  many  inverse  problems  in  geophysics  both  of  fundamental  interest  such  as  the 
reconstruction  of  the  Earth  structure  from  knowledge  of  elastic  ,  gravimetric  and  geomagnetic 
data  or  of  more  applied  nature  such  as  the  use  of  similar  data  in  geophysical  prospecting  to 
localize  gas  or  oil. 

In  this  paper  we  consider  an  inverse  problem  in  gravimetry.  Let  x,  y,  z  be  cartesian  coordinate 
in  the  three  dimensional  space,  given  a  mass  density  p(x,y,z)  in  a  certain  region  fi,  the  direct 
problem  of  gravimetry  consists  in  finding  the  gravitational  potential  V(x,y,z)  generated  by  p, 
that  is  in  solving  the  Poisson’s  equation: 

AK  = -4* fp  in  n  (1.1) 

where  A  =  +  is  the  laplacian,  /  =  6.67- 10~ 8 cm3  • g/atc 2  is  the  gravitational  constant 

and  suitable  boundary  conditions  on  the  boundary  of  fl,  dfl,  are  specified  Vice  versa  given  the 
gravitational  potential  V(z,y,z)  in  a  certain  region  Q,  the  inverse  gravimetry  problem  consists  in 
finding  the  mass  density  p(x,y,z)  that  generates  V,  that  is  in  solving  the  inverse  Poisson’s  equation 

A-1p=-^V  infl  (1.2) 

where  A-1  is  the  inverse  laplacian  with  the  appropriate  boundary  conditions  on  <90. 

Almost  all  inverse  problems  and  in  particular  the  inverse  gravimetry  problem  are  ill  pcsed, 
that  is,  to  an  arbitrarely  small  perturbation  of  the  data  V  can  correspond  an  arbitrarely  large 
perturbation  of  the  solution  p. 

The  great  difficult  cf  solving  the  inverse  gravimetry  problem  in  practical  situations  is  due 
to  its  ill  posedness  since  the  values  of  the  gravitational  field,  that  is  the  data,  are  obtained  from 
experiments  and  so  are  affected  by  errors.  In  order  to  avoid  this  difficulty  and  restore  a  well 
behaved  dependence  of  the  solution  p  from  the  data  V  it  is  very  useful  to  introduce  some  a  priori 
constraints  that  o  should  satisfy. 

In  [3],  [4)  [5]  using  these  ideas  the  inverse  gravimetry  problem  has  been  reduced  to  a  linear 
programming  problem  in  a  way  that  we  will  see  later.  In  this  paper  using  this  linear  programming 
formulation  of  the  inverse  gravimetry  problem  we  study  the  northern  San  Francisco  craton  granite 
in  Brazil  on  the  basis  of  measurements  of  the  residual  Bouguer  anomaly  of  the  gravitational  field 
taken  from  [6]  and  [7],  see  Fig.l.  From  these  data  we  reconstruct  a  two  dimensional  section  of  the 
mass  density  along  the  BB’  segment  of  Fig.l.  This  section  is  about  210Km  long  and  20 Km  deep, 
see  Fig. 3.  Moreover  we  reconstruct  a  three  dimensional  section  of  the  mass  density  in  th’  region 
bounded  by  the  dotted  rectangle  of  Fig.l.  This  section  on  the  surface  of  the  Earth  is  a  rectangle 
of  about  240fifm  x  189Km  and  is  about  20 Km  deep,  see  Fig. 4. 

We  have  found  that  in  the  center  of  this  three  dimensional  .region  there  is  a  body  of  granite 
of  mass  density  p  =  2.57 g/cm3  while  the  average  mass  density  of  this  region  is  p  =  2.67 j/cm3. 
This  body  in  the  BB'  segment  is  surfacing  for  about  30 Km  and  is  about  67  Km  long  on  its  bottom 
which  is  about  16 Km  deep,  see  Fig. 3.  These  results  are  confirmed  by  the  three  dimensional 
reconstruction  of  p  (Fig. 4)  and  they  are  in  good  agreement  with  the  results  obtained  by  I'ssami 
and  Bott  in  [6],  The  granite  body  found  by  Ussami  and  Bott  in  [6]  is  a  little  bit  smaller  than  ours. 
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In  section  2  we  describe  our  mathematical  model  and  we  reduce  the  inverse  gravimetry  problem 
to  a  linear  programming  problem.  In  section  3  we  present  the  data  that  we  have  used  and  the 
results  obtained  in  the  northern  San  Francisco  craton  granite  (Brazil)  using  the  mathematical 
model  described  in  section  2. 

2 .  The  Mathematical  model 

Here  we  will  follow  the  work  of  Safon,  Vasseur  and  Cuer  [4],  it  is  well  known  that,  using  the 
Green’s  function  G(x,y,z )  [8],  the  equation  (1,2)  can  be  reformulated  as  an  integral  equation  as 
follows 

/  [  0{xtytztxtty9tzt)p(xttytiz,)dxtdyldzt  =  V(xiy}z)  (2.1) 

J  n 

Let  us  cut  the  domain  fl  in  a  great  number  N  of  cubic  domains  uni  the  function  p(x,  y,  z)  will  be 
approximated  on  each  cube  <u„  by  a  constant  pn  and  the  Green’s  function  G  will  be  approximated  by 
a  matrix  £  =  ((j,y))  ;  i  =  1,2,  ••  • ,  M  ;  j  —  1,2,  •  •  • ,  N .  Finally  let  V)  ,  i  =  1,2,  ■  •  • ,  M  (M  <  N) 
be  the  measurements  of  the  gravitational  field  V(x,y,  z)  in  x, ,  y, ,  z,  that  belongs  to  the  cubic 
domain  w,,«  =  1,2,---,  M.  So  that  equation  (2.1)  can  be  discretized  as  follows: 

£*>P>=V'  ,  •  =  l,2,---,M  (2.2) 

j=i 

where  g,,  =  /  J  C(x,,  y, ,  Zi,  z' ,  y',  z')dx'dy'dz'.  In  general  is  natural  to  assume  that  the  data  V, 
are  available  only  on  cubic  domains  w,  that  have  one  side  on  the  accessible  surface  of  the  region 
fl,  so  that  M  <  N ,  see  Fig. 4  and  the  linear  system  (2.2)  is  underdetermined.  Moreover  since  the 
data  V,  ,  i  =  1,2,  ■  ■  ■  ,  M  are  known  only  with  a  certain  error  t,  >  0  ,  »  =  1,2,  •  •  • ,  M  we  substitute 
the  system  of  linear  equation'  (2.2)  with  the  system  of  linear  inequalities 

,v 

V.-Ci  <  S  V.+e,  ;  t  =  1,2,  -  ■  • ,  Af  (2.3) 

y=i 

The  mass  density  p  is  greater  or  equal  to  zero,  so  that  we  may  impose 

p,  >0,  j=  1,2,  (2.4) 

or  more  realistically  in  a  geophysical  problem 

pT"  <  P,  <  P™  ,  3  =  IX--,  If  (2.5) 

where  p™"  and  ,  j  =  1,2, •••,1V  are  the  minimum  and  the  maximum  values  for  p}  assigned 
on  the  basis  of  physical  intuition.  Both  the  linear  system  (2.2)  or  the  more  realistic  system  of 
inequalities  (2.3),  (2.5)  if  compatible  have  in  general  an  infinite  number  of  solutions. 

In  order  to  restore  uniqueness  we  introduce  a  functional  to  be  optimized  on  the  set  of  points 
satisfying  (2.3),  (2.5) ,  that  is 

L 

j  =  YLpk’  v°l(“k.) 

•=i 


(2.6) 
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where  voJ(vk,)  is  the  volume  of  the  cube  ,  so  that  J  is  the  total  mass  of  the  body  0  or  of  some 
part  of  it. 

The  linear  programming  problem  that  we  consider  is  the  following  one:  minimize  the  functional 
J  given  by  (2.6)  subject  to  the  constraints  (2.3),  (2.5). 

In  section  3  we  will  also  use  a  two-dimensional  model,  in  this  case  we  assume  that  the  region 
fl  has  as  infinite  extension  in  the  horizontal  x  coordinate  and  that  p  is  a  function  of  z,  the  vertical 
coordinate,  and  y  only.  Reasoning  as  in  the  three-dimensional  case  with  an  appropriate  Green’s 
function  G  and  data  V,  ,  i  =  l,2,---,Af  it  is  easy  to  obtain  a  linear  programming  problem 
analogous  to  the  one  obtained  for  the  three-dim'nsional  model  given  by  (2,3),  (2.5),  (2  6). 


3.  Results  and  conclnsions 

Here  we  will  consider  two  special  applications  of  the  method  described  in  section  2.  That  is  a 
two  dimensional  vertical  section  and  a  three  dimensional  slice  of  a  large  region  in  the  state  of  Bahia 
in  Brazil.  In  our  computations  the  gravitational  anomaly  is  measured  in  mg  ala  and  the  distances 
in  Km.  The  gravitational  anomalies,  that  is  the  residual  Bouguer  anomalies,  are  reported  in  Table 
1  for  the  three-dimensional  problem  and  in  Fig. 2  for  the  two-dimensional  problem.  Moreover  we 
assume  that  the  average  mass  density  in  the  regions  considered  is  po  —  2.67 g/ cm3  The  data  of  table 
1  and  Fig. 2  have  been  regularized  in  order  to  avoid  numerical  instability  due  to  the  ill  posedness 
of  the  problem.  In  particular  the  data  of  Fig. 2,  that  are  taken  along  the  BB’  segment  of  Fig  1, 
are  highly  irregular  outside  the  central  region  of  Fig. 2  where  the  maximum  of  the  gravitational 
anomaly  is  attained  so  that,  as  shown  in  Fig. 2,  the  data  have  been  regularized  outside  this  region. 
For  the  three-dimensional  case  the  numbers  in  bold  face  of  Table  1  are  real  gravitational  anomalies, 
these  are  the  data  taken  in  the  rectangle  AA’C’C  of  Fig.l,  while  the  numbers  not  in  bold  face  in 
Table  1,  that  is  the  data  of  the  region  between  the  dotted  rectangle  and  AA’C’C  of  Fig.l,  have 
been  regularized  solving  some  suitable  direct  problem.  In  order  to  choose  the  functional  J  given 
by  (2.6)  we  have  used  an  empiric  formula  [9j  that  estimates  the  maximum  deepness  Az  of  a  body 
that  generates  a  certain  gravitational  anomaly  V,  ,  i  =  1,2,  •  •  • ,  M  on  the  surface,  that  is 


Az  =  0.65 


max,  IV,  | 
max,  |  v  V. 


(3.1) 


where  vK  is  a  numerical  approximation  to  the  gradient  of  V,  expressed  in  mgala/ Km.  With  our 
data  we  obtain  A z  ~  16 Km.  So  that  on  the  basis  of  this  estimate  for  A z  and  of  the  conclusion  of 
[6]  that  the  body  is  surfacing  we  have  chosen 

J  =  ^2  Pi  voluj  (3  2) 

{l  (  <3Km.<t<l6Km.) 

In  particular  in  the  two-dimensional  case  along  the  BB’  segment  of  Fig.l  we  have  considered  37 
data  points  taken  on  the  surface  at  a  distance  of  7.5 Km  one  from  the  other.  Moreover  we  assume 
that  the  body  we  are  looking  for  is  localized  in  the  region  z  <  2 OKm.  So  that  the  rectangle  we 
are  working  on,  as  shown  in  Fig. 3,  is  (7.5 Km  x  37  =  277.5 Km)  x  20 Km  and  has  been  divided  in 
small  rectangles  w,  ,  j  =  1,2, •••,370  where  u,  is  7.5/fm  x  2 Km  (see  Fig  3)  and  j  runs  from  left 
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to  right  and  from  top  to  bottom  from  1  to  370.  The  bound  (2.5)  assumed  on  p,  ,  j  =  1, 2,  •  •  • ,  370 
have  been  chosen  so  that  the  feasible  region  is  as  small  as  possible. 

The  gravimetric  problem  that  has  been  solved  for  the  two  dimensional  case  is  the  following 
one:  find  the  values  of  py  ,  j  =  1, 2,  •  •  • ,  370  that  minimise 


2«e 

J  =  ^  p'VoltjJ , 

J=1 


(3.3) 


subject  to 


370 

Vi  p>  - Vx 

2.57 g/cm3  <  p,  <  2.69 g/cm3 
2.57 g/cm3  <  p,  <  2.675 g/cm3 
2.57 g/cm3  <  p,  <  2.67g/cm3 


i  =  l, 2, •••,37 

;  l  <  j  <  37 
;  38  <  j  <  74 
i  75  <  J  <  370 


(3.4) 


where  V)  =  4  +  7.  with  5,  =  2.67g/cm3  £*I°i  p,y  and  7,  the  anomalies  shown  in  Fig.2  while 
t,  =  O.O57,.  This  linear  programming  problem  (3. 3), (3.4)  has  been  solved  using  the  FORTRAN 
package  [10]  and  the  solution  found  is  shown  in  Fig. 3. 

In  the  three  dimensional  problem  the  dotted  rectangle  of  Fig.  1  is  of  dimension  240 Am  in  the 
BB’  direction  (y  direction)  times  189Am  in  the  AC  direction  (z  direction).  This  dotted  rectangle 
has  been  divided  in  16  X  7  =  112  small  rectangles  of  dimension  15  Am  in  the  BB’  direction  times 
27 Km  in  the  AC  direction,  the  gravitational  anomaly  data  7,  are  taken  in  the  central  points  of 
these  small  rectangles  {xit  y3 ,0)  ,  i  =  1,2,  •  •  •  ,7  ,  j  —  1,2, •••,16. 
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The  parallelepiped  considered  is  the  one  of  base  the  dotted  rectangle  of  Fig.l  and  deepness  20Km. 
This  parallelepiped  has  been  divided  in  7  x  16  x  5  =  560  small  parallelepipeds  coy  ,  j  =  1,2,  •  •  • ,  560 
of  dimensions  27 Km  x  15Km  x  4fCm  in  the  z,y,z  directions  respectively  (see  Fig.  4)  and  j  runs 
as  shown  in  Fig.  4. 

The  gravimetric  problem  that  has  been  solved  with  the  FORTRAN  package  [10]  is  the  following 
one:  find  the  values  of  p}  ,  j  =  1,2, •••,560  that  minimize  J  given  by  (3.2)  subject  to: 
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560 

v.  -  £.  <  <  V.  +  e.  .•  =  1.2.  — ,  112 

j=i 

2.57g/cms  <  p;  <  2.697 j/cm!  ;  1  <  jf  <  112 

2.57 g/em3  <  p,  <  2.672g/cm3  ;  113  <j<  224 
2.57 g/em3  <  p,  <  2.61g/cm 3  ;  225  <  j  <  560 

where  e,  and  V,  =  5,  +  7,  are  defined  aa  in  the  two-dimensional  case.  The  results  obtained  are 
shown  in  Fig.4. 

As  shown  in  Fig.3  for  the  two-dimensional  problem  we  have  found  a  granite  body  about  16 Km 
deep  with  density  2.57 j/cm3,  that  is  0.1  g/cm3  smaller  that  the  surrounding  medium,  this  body 
is  about  30 Km  long  on  the  surface  and  67 Km  long  on  the  bottom.  The  results  of  the  three- 
dimensional  study,  shown  in  Fig.4,  confirm  the  results  of  Fig.3  This  results  are  in  qualitative 
agreement  with  the  ones  obtained  by  Ussami  and  Bot  in  [6j. 
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FIGURE  CAPTIONS 


Fig.l  A  large  region  in  the  state  of  Bahia  (Brazil)  is  considered  between  40-44  degrees  of 
longitude  and  10-15  degrees  of  latitude.  The  iso-Bouguer  anomaly  lines  are  shown.  The 
BB’  direction  is  the  one  considered  in  the  two-dimensional  problem.  The  dotted  rectangle 
is  the  region  considered  in  the  three-dimensional  problem.  In  the  three-dimensional 
problem  measured  data  have  been  used  in  the  AA’C’C  rectangle,  regularized  data  have 
been  used  in  the  remaining  region. 

Fig.2  The  measured  residual  Bouguer  anomaly  in  the  BB’  segment  and  the  regularized  one  is 
shown. 

Fig.3  The  mass  density  distribution  obtained  from  the  data  of  Fig.2  solving  the  linear 
programming  problem  (3. 3), (3. 4)  in  a  vertical  section,  20 Km  deep, along  the  BB’ 
direction  is  shown. 

Fig. 4  The  mass  density  distribution  obtained  from  the  data  of  Table  1  solving  the  linear 
programming  problem  (3. 2), (3. 5)  is  shown. 
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Abstract.  Serial  and  parallel  successive  overrelaxation  (SOR)  methods  are  proposed  for 
the  solution  of  the  augmented  Lagrangian  formulation  of  the  dual  of  a  linear  program.  With 
the  proposed  serial  version  of  the  method  we  have  solved  linear  programs  with  as  many  as 
125,000  constraints  and  500,000  variables  in  less  thac  72  hours  on  a  MicroVax  II.  A  parallel 
implementation  of  the  method  was  carried  out  on  a  Sequent  Balance  21000  multiprocessor 
with  speedup  efficiency  of  over  65%  for  problem  sizes  of  up  to  10,000  constraints,  40,000 
variables  and  1,400,000  nonzero  matrix  elements. 
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1.  Introduction 

In  [8,  9,  10,  12,]  successive  overrelation  methods  are  proposed  for  solving  the  dual  of 
the  problem  of  finding  the  least  2-norm  solution  of  a  linear  program.  This  leads  to  an 
exterior  penalty  formulation  of  the  dual  of  the  original  linear  program  with  the  interesting 
property  that  the  penalty  parameter  need  not  approach  zero  in  order  to  obtain  an  <*xact 
solution  of  the  primal  linear  program  [1,  14].  Thus  the  penalty  parameter  need  only  be 
less  than  a  certain  threshold  value  in  order  to  obtain  an  exact  solution  to  the  primal  linear 
program.  However,  the  penalty  parameter  must  approach  zero  in  order  to  obtain  a  solution 
to  the  dual  problem.  Although  this  approach  has  been  used  effectively  in  conjunction 
with  successive  ov^rrelaxation  methods  both  on  serial  [10]  and  parallel  machines  [11,  12], 
we  propose  here  the  use  of  an  augmented  Lagrangian  on  the  dual  problem  instead  of 
an  exterior  penalty  function  in  order  to  alleviate  the  twin  difficulties  of  determining  the 
threshold  value  of  the  penalty  parameter  required  for  an  exact  primal  solution,  and  the 
need  for  the  penalty  parameter  to  approach  zero  in  order  to  obtain  a  dual  solution.  The 
first  proposal  for  using  an  augmented  Lagrangian  formulation  for  solving  linear  programs 
was  made  in  [22].  In  [18]  Polyak  and  Tretiyakov  made  the  remarkable  discovery  that  after 
a  finite  number  of  steps  of  the  augmented  Lagrangian  algorithm,  an  exact  solution  to 
the  primal  and  dual  linear  programs  is  obtained.  In  [3]  Golshtein  proposed  a  projected 
Gauss-Seidel  method  in  conjunction  with  an  augmented  Lagrangian  formulation  and  gave 
computational  results  for  linear  programs  with  sizes  up  to  352  variables  and  166  constraints. 
No  convergence  proofs  of  the  projected  Gauss-Seidel  method  wcs  given  in  [3],  nor  of  the 
closely  related  iterative  method  of  Syrov  and  Churkreidze  in  [IS],  We  propose  here  the 
use  of  a  projected  successive  overrelaxation  method  in  conjunction  with  an  augmented 
Lagrangian  formulation.  The  convergence  of  the  projected  SOR  method  established  in 
[7]  is  general  enough  to  cover  both  a  serial  and  a  parallel  implementation  of  the  method. 
Since  SOP.  methods  are  inherently  serial  in  n;  tura,  their  parallelization  is  not  a  routine 
matter.  In  [11,  13]  two  related  methods  were  proposed  for  the  parallelization  of  SOR 
methods.  The  more  recent  method  [13]  utilizes  an  unreduced  relaxation  factor  interval  of 
(0,  2)  which  we  shall  employ  here  with  an  augmented  Lagrangian  algorithm  for  the  dual 
linear  program. 

The  paper  is  organized  as  follows.  In  Section  2  we  give  the  necessary  theoretical 


background  and  convergence  results  for  the  proposed  augmented  Lagrangian  method  ap¬ 
plied  to  the  dual  linear  program.  Tn  Section  3  we  describe  a  serial  SOR  implementation 
of  the  method  and  establish  its  convergence.  In  Section  4  we  describe  our  parallel  SOR 
implementation  and  in  Section  5  we  present  computational  results  for  both  the  serial  and 
parallel  methods. 

We  briefly  describe  our  notation  now.  For  a  vector  x  in  the  n-dimensional  real  space 

Rn ,  x+  will  denote  the  vector  in  Rn  with  components  (x+).  =  max  {x,,  0 } ,  i  =  1 . n. 

The  scalar  product  of  two  vectors  x  and  y  in  Rn  will  be  simply  denoted  by  xy .  For 

n 

1  <  p  <  oo ,  the  p-norm  ^  |x,  j^)  of  a  vector  in  R"  will  be  denoted  by  ||  x ! j p  .  For  the 

i=t 

2-norm  the  subscript  2  will  be  dropped.  R"  will  denote  the  nonnegative  orthant  or  the  set 
of  points  in  Rn  with  nonnegative  components,  while  Rmxn  will  denote  the  set  of  all  m  x  n 
real  matrices.  For  .4  £  Rmxn,  AT  will  denote  the  transpose.  .4,  will  denote  the  ith  row, 

A,j  the  element  in  row  i  and  column  j ,  and  for  I  C  {1 . m},  J  C  {l . n},  .4;  will 

denote  the  submatrix  of  .4  with  rows  .4,,  t  £  /  ,  while  Aij  will  denote  the  submatrix  of  .4 
with  elements  ,4;;-,  i  £  I,  j  £  J.  Similarly  for  x  ~  Rn  and  I(  C  {1, .  . .  ,  n},  x/(  will  denote 

x,,  i  £  Ie.  The  set  {/t)  /2, . . . ,  IK}  is  said  to  be  a  consecutive  partition  of  {1 . n}  if 

it  is  a  partition  of  {1, - n}  such  that  i  <  j  for  i  £  j  e  It+\  and  (  =  1 _ ,  Jfc  -  1 . 

Here  and  throughout  the  symbols  :  =  and  =:  denote  definition  of  the  term  on  the  left  and 
right  sides  of  each  symbol  respectively. 
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2.  Theoretical  Background 

We  consider  the  linear  program 

(2.1)  min  ex  subject  to  Ax  >  b,  x  >  0 
where  c  €  Rn ,  6  6  Rm  and  A  £  i?mxn  and  its  dual 

(2.2)  max  bu  subject  to  v  =  —ATu  +c  >  0 

(u,u)>0 


For  simplicity  we  exclude  trivial  constraints  with  Aj  =  0.  In  [8,  9]  the  exterior  penalty 
problem  associated  with  the  dual  problem  (2.2) 

(2.3)  max  eiu- i||ATu+v-c|]2 

(u,v)>  0  6 


was  solved  by  an  SOR  procedure  for  a  sufficiently  small  value  of  the  penalty  parameter 
e  to  obtain  (u(e),  v(e)) .  The  unique  least  2-norm  solution  x  of  the  linear  program  (2.1) 
was  obtained  by  using  the  equation 

(2.4)  x(f )  =  ^  (Aru(e) -f  u(e)  -  c) 

which  relates  an  optimal  solution  (u(c),  v(e))  of  the  dual  penalty  problem  (2.3)  and  the 
unique  solution  x(e)  of  the  corresponding  quadratic  primal  problem  [5] 

(2.5)  min  cx  +  ~  xx  subject  to  Ax  >  b,  x  >  0 

In  particular  it  follows  [14]  that  the  least  2-norm  solution  x  of  the  linear  program  (2.1)  is 
related  to  x(e)  of  (2.4)  by 


(2.6) 


x  =  x(e)  for  £  6  (0,  e]  for  some  c  >  0 


Thus  the  penalty  parameter  e  of  the  dual  penalty  problem  (2.3)  is  the  perturbation  param¬ 
eter  of  the  perturbed  primal  problem  (2.5).  In  order  to  avoid  possible  difficulties  associated 
with  determining  the  threshold  value  £  we  consider  instead  of  the  exterior  penalty  problem 
(2.3)  the  augmented  Lagrangian  associated  with  the  dual  linear  program  (2.2) 

L(u,  v,  x,  7)  :=  bu  —  [[ATu  4-  v  —  c||‘  —  x(ATu  +  v  —  c) 


(2.7) 


26 


It  is  a  standard  result  [6,  2,  19]  that  for  any  7  >  0,  a  primal-dual  solution  (i,u,v)  of 
(2.1)-(2.2)  is  equivalent  to  a  stationary  point  of  the  following  saddlepoint  problem  of  (2.7): 
Find  an  (z,-u,u)  £  Rn  x  R™  x  such  that  for  all  x  £  Rn  and  all  (u,u)  £  R+  x  R™, 

(2-8)  L(u,v,  x,  7)  <  L{i 1,0, 1,7)  <  L(u,v,x,  7) 

The  standard  augmented  Lagrangian  algorithm  [2,  19]  consists  of  a  maximization  step  in 
the  (u,v)  space  R™  x  R+  followed  by  an  unconstrained  gradient  descent  step  in  the  x 
space  Rn .  In  particular  we  have  the  following. 


2.1  Augmented  Lagrangian  Algorithm 

Start  with  any  x°  £  Rn .  Having  x1  determine  x‘+l  as  follows 


(2.9) 


(a)  L(u\  v',  x‘,  71)  =  max  L(u,  v,  x‘,  7') 

(u,u)>0 

(b)  x,+1  =  x1'  -  i-V^u Sv'.xSy)  =  X1  +  i(Aru'  +  v'  -  c) 

7  7* 


where  {7*}  is  any  bounded  sequence  of  positive  numbers. 

For  this  iterative  linear  programming  algorithm  Polyak  and  Tretiyakov  [IS]  have  given 
the  following  important  finite  convergence  theorem. 


2.2  Augmented  Lagrangian  Algorithm  Finite  Termination  Theorem  [18] 

For  any  bounded  positive  sequence  {71}  and  x°  £  Rn,  Algorithm  2.1  is  finite,  that  is 
there  exists  an  integer  k  such  that  (xk,vk,vk)  solve  the  dual  linear  programs  (2.1  )-(2.2). 
Furthermore  for  each  x°  there  exists  a  7  >  0  such  that  for  0  <  70  <  7,  the  method  will 
terminate  in  one  step,  that  is  (x’.u1,!;1)  will  solve  the  dual  linear  programs  (2.1)-(2.2). 
Note  that  in  the  above  theorem,  two  exact  maximizations  over  the  (u,u)  space  R™  x 
are  required  in  order  to  obtain  (x,+1,u,+  1,u*+l)  from  x'. 

We  note  that  since  by  duality  theory  [oj 


(2.1°)  (maxa  L(u,  v,  x,  7)  =  mm{cx+  1 1|2  -  x||2j.4x  >  6,  x  >  0}  =:  <p(x) 

the  Augmented  Lagrangian  Algorithm  2.1  is  equivalent  to  the  following  gradient  method 
applied  to  the  proximal  point  function  <b(x) 

f  xi+1  =  *•'  -  i- Vd(x')  =  Prox(xi) 


(2.11) 


where  Prox(x')  is  the  solution  of 
^  mm  {  cz  +  2-.\\z~  x'  || 2  j  Ax  >  b,  z  >  0} 
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Bertsekas  [l]  and  Rockafellar  [21]  also  obtain  finite  termination  for  (2.11)  from  proximal 
point  theory  considerations  for  the  linear  programming  case. 

With  this  background  we  are  prepared  now  to  state  and  prove  the  convergence  of  our 
serial  and  parallel  SOR  algorithms. 
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3.  Serial  Successive  Overrelaxation  Algorithm 

The  proposed  serial  algorithm  consists  of  applying  the  projected  SOR  method  of  f7|  to 
the  maximization  step  (2.9a)  of  the  Augmented  Lagrangian  Algorithm  2.1  for  a  decreasing 
sequence  of  positive  numbers  {71}-  It  follows  from  the  Finite  Termination  Theorem  2.2 
and  a  theorem  of  Pang  (16,  Theorem  3.1],  that  for  any  x'  the  projected  SOR  method 
will  generate  a  sequence  of  points  converging  to  an  x,+  1  that  solves  the  primal  linear 
program  (2.1),  provided  that  Y  is  sufficiently  small.  There  are  no  easily  implementable 
and  theoretically  justifiable  ways  of  determining  how  to  choose  71  sufficiently  small  [20,  4] , 
however  we  shall  prescribe  some  computationally  effective  ways  for  doing  that. 

The  serial  projected  SOR  method  has  been  proposed  (7,  9]  for  solving  the  quadratic 
minimization  problem 


(3.1) 


min  9(z )  :=  min  -  zMz  +  qz 

z> o  t> o  2 


where  M  C  Rkxk  is  symmetric  and  positive  semidefinite.  This  is  precisely  problem  (2.9a) 
of  Algorithm  2.1  if  we  make  the  identifications 


(3  2) 


M:  =  —  | 

A 

U- 

1 

1 

o- 

,  Z  \  — 

u 

Y 

I 

L  x.  -  4  J 

V 

x' - 

7’ 


The  projected  SOR  algorithm  consists  of  the-following. 

(33)  4*'  -  (*:  . . 4)) 

W  €  (0,  2),  j  =  1,.  ..,k 

More  specifically  for  9(z)  =  i  zMz  +  qz  we  have  the  following. 

3.1  Serial  SOR  Algorithm  for  min  i-zMz  +  qz 

*>o  2 

Choose  2°  6  iZ*  ,  w  €  (0,  2).  Having  z‘  compute  z<+1  as  follows 


+ 


(3.4) 


3-1 


r  =  (*>-  (  e  ^+i + e  + «) 


t= i 
for  ;>1 


1=3 


We  are  ready  to  state  and  establish  the  convergence  of  our  augmented  Lagrangian 
serial  SOR  algorithm. 

3.2  Augmented  Lagrangian  Serial  SOR  Algorithm 

Let  {t1}  I  7  for  some  y  >  0,  let  {<f‘}  j  0  and  let  x°  £  R" .  Having  x‘  determine 
x,+  1  as  follows: 

(a)  Apply  the  Serial  SOR  Algorithm  3.1  to  solve  (2.9a)  with  the  identifications  (3.2),  and 
let  (u‘(t),  u‘(t))  be  the  t  iterate  of  this  SOR  algorithm.  Stop  if  for  some  t  =  t'  the 
following  inequality  is  satisfied. 

|u‘(f)V uX(u'(t),  x',  7')|  +  It ;‘(t)V„L(u1(t),  v‘{t),  x‘,  7’)  I 

(3  5) 

+  j| (VaL(u‘(t),  id(t),  x\  7')  +  J  +  |j(v„L(u'(<),  v‘(t),  x',  7'))  J  <  S' 

(b)  Set  x,+I  =  x'(f')  where 

(3.6)  x‘(f)  :=  x-  +  ~(ATu'(t)  +  v'(t)  -  c) 

7* 


3.3  Augmented  Lagrangian  Serial  SOR  Convergence  Theorem 

Let  {71}  l  i  >  0  be  a  sufficiently  rapidly  decreasing  sequence  of  positive  numbers 
and  7  sufficiently  small.  Then  either 

(a)  For  some  integer  k,  the  sequence  (x*(<)}  converges  to  an  x*  that  solves  the  linear 
program  (2.1),  or 

(b)  For  each  subsequence  of  {(x‘,  u'(t'),  u'(f'))}  converging  to  some  (x,u,u),  the  cor¬ 
responding  subsequence  {x,+1  =  x'(t')}  converges  to  an  x  such  that  i  solves  the 
linear  program  (2.1).  If  x  =  x,  then  (u,v)  solves  the  dual  linear  program  (2.2). 

Proof  Either  the  inequality  (3.5)  of  the  algorithm  is  satisfied  at  each  iteration  i  for  some 
t 1  or  not.  Accordingly  we  have  the  alternatives  (b)  and  (a)  below  respectively. 

(a)  For  some  iteration  i  =  k  the  inequality  (3.5)  is  never  satisfied.  Hence  by  Pang's 
Theorem  3.1  [16],  since  L(u,v,xk  ,yk)  is  by  duality  theory  bounded  above  for  (u,u)  > 
0  by  cx  +  7*||x  —  x*||2  for  any  x  >  0  such  that  Ax  >  b,  it  follows  that  the  sequence 


{**«)}  - 


xk  +  ATuk(t)  +  vk{t)  -  c 


},  <  =  0,1,2,. 


where  t  is  the  SOR  iteration  index,  converges  to  a  vector  x*  defined  by 
(3.3)  +  xj: 

for  some  ( uk,vk )  which  solves  max  £(u,  u,  xk ,  ik).  Note  however  that  {uk(t)%  vk(t)} 
need  not  converge  to  (uk,vk).  If  7*  is  sufficiently  small  (and  this  is  what  is  meant  by 
requiring  that  {7'}  decreases  sufficiently  rapidly)  it  follows  by  Theorem  2.2  that  xk+l  — 
xk,  is  a  solution  of  the  primal  linear  program  (2.1).  (Note  that  [uk,vk)  need  not  be  a 
solution  of  the  dual  linear  program  (2.2).  To  obtain  such  a  dual  optimal  (u*,  uk )  we  need 
to  solve  min  £(  (u,  v,  x k,  yk)  exactly.  See  Corollary  3.5  below.) 

(u,v)>0 

(b)  If  inequality  (3.5)  holds  for  each  t  for  some  t\  then  since  { 5* }  j.  0  and  {7'}  J, 

7  >  0,  we  have  that  for  any  subsequence  of  {(x',u‘(t*),  u'(t‘))}  converging  to  some 

(z,u,v),  the  point  (u,u)  solves  the  problem  min  £(u,u,x,  7)  (because  (<5‘ }  J.  0), 

(u*v)>0  ' 

and  moreover  the  corresponding  subsequence  |x,+1  }  defined  by  (3.6)  and  part  (b)  of 
Algorithm  3.2  converges  to  an  i  defined  by 


x  x  + 


A1  u  +  v  —  c 


If  7  is  sifficiently  small,  then  it  follows  by  Theorem  2.2  that  x  solves  the  linear  program 
(2.1).  If  in  addition  x  =  x,  then  (u,v)  is  feasible  for  the  dual  Linear  program  (2.2)  and  is 
also  optimal  because  bu  =  cx.  I 

It  is  useful  to  point  out  that  when  the  sequence  {(u'(t),  v‘(t)) }  of  Algorithm  3.2  is 
bounded  then  it  has  an  accumulation  point  and  by  Theorem  2.1  [7],  each  such  accumulation 
point  satisfies  the  optimality  conditions  for  (2.9a).  Consequently  for  each  i  inequality  (3.5) 
of  Algorithm  3.2  is  satisfied  after  a  finite  number  of  steps  of  part  (a)  of  the  algorithm. 
However  by  Theorem  2(ii)  of  [9]  we  have  that  if  the  linear  program  (2.1)  satisfies  the  Slater 
constraint  qualification,  then  the  sequence  {(u‘(t),  u'(t))}  is  indeed  bounded.  Therefore 
we  have  the  following. 

3.4  Augmented  Lagrangian  Serial  SOR  Convergence  Corollary 

Let  the  linear  program  (2.1)  satisfy  a  Slater  constraint  qualification,  that  is  Ax  >  b 
for  some  x  >  0.  Then  Theorem  3.3  holds  with  outcome  (b). 
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Another  useful  observation  follows  from  the  fact  [18,  Lemma  1],  [19]  that  minimizing 
the  augmented  Lagrangian  £(u,u,x,  7)  with  an  optimal  value  of  r  and  any  7  >  0  gives  a 
solution  to  the  dual  linear  program  (2.2).  Hence  we  have  the  following. 


3.5  Dual  LP  SOR  Solution  Corollary 

Under  the  assumptions  of  Theorem  3.3  a  solution  to  the  dual  linear  program  (2.2)  can 
be  obtained  for  either  outcome  (a)  or  (b)  of  Theorem  3.3  by  solving  respectively 

(a)  min  i(u,u,x*,7l) 

(u,u)>v 

or 

(b)  min  £(u,u,z,7) 

(u,v)>¥' 

We  note  immediately  that  we  have  left  open  the  procedure  by  which  the  sequence 
{7'}  is  decreased.  This  is  an  inherent  theoretical  difficulty  that  arises  when  using  inexact 
minimization  in  the  subproblems  of  proximal  point  or  augmented  Lagrangian  algorithms. 
Thus  the  approximate  minimization  criteria  of  [21,  Criteria  A,  B,  A' ,  B'}  are  not  imple- 
mentable  for  our  problem,  while  the  assumptions  of  [2,  Section  2.5]  are  not  verifiable  for  our 
problem.  Computationally  we  have  overcome  this  difficulty  by  using  the  following  scheme, 
often  used  for  updating  the  penalty  parameter  in  augmented  Lagrangian  algorithm 


(3.10) 


7'  if  ||x‘+1  -x'||  <^||x'  -x,-‘||,  0<  »  <  1 
uf‘  otherwise,  0  <  u  <  1 


This  scheme  works  effectively  for  the  solution  of  very  large  sparse  linear  programs  as  our 
computational  results  indicate. 
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4.  Parallel  Successive  Overrelaxation  Algorithm 

The  key  to  our  parallel  SOR  algorithm  is  the  use  of  the  parallel  gradient  projection 
SOFt  (GP-SOR)  that  we  proposed  in  (13)  for  the  solution  of  (3.1)  and  which  we  outline 
below.  Partition  the  matrix  M  of  (3.1)  into  r  contiguous  horizontal  blocks  as  follows: 

'A/,, ' 

Mi, 

(4.1)  M  =■ 

-Mi.. 

where  the  blocks  Mr,  correspond  to  the  variables  ,  and  (Jj,  Ii, . . . ,  Ir}  is  a  consecutive 
partition  of  {1, 2, . . . ,  k)  .  Now  partition  Mi,  as  follows 

(4.2)  Mi,  =:  (A//,  Mj.f.] 

where  Ij  is  the  complement  of  I}  in  {1, 2, . . . ,  k} .  Thus  Mi,i,  is  a  principal  square 
submatrix  of  M  with  elements  Mat,  s  €  Ij  and  t  6  Ij.  We  further  partition  Mi, i,  into 
the  sum  of  its  strictly  lower  triangular  part  ,  its  diagonal  part  D^i,  and  its  strictly 
upper  triangle  part  Ufji,  as  follows 

(4.3)  Mi, I,  =:  Li,!,  +  L>r,i,  4-  U,,i, 


Now  define  a  block  diagonal  matrix  K  as  follows 


(4.4) 


K:  = 


Li, I, 


Li. i.  J 


where  each  L[,i,  is  the  strictly  lower  triangular  part  of  M[,i ,  .  An  SOR  algorithm  can 
now  be  performed  for  each  row  block  Ij,  j  =  1, . . .  ,r,  simultaneously,  that  is  in  parallel. 
Note  that  the  block  diagonal  matrix  K  replaces  the  traditional  strictly  lower  triangular 
matrix  of  the  serial  SOR.  Specifically  we  have  the  following. 

4.1  Parallel  GP-SOR  Algorithm  for  (3.1) 

Let  {Ii,  Ii, . . . ,  Ir}  be  a  consecutive  partition  of  {1,2,... ,fc},  let  £  be  a  positive 
diagonal  matrix  in  and  let  z°  >  0.  For  t  =  0, 1,  2, ... ,  do  the  following 


I 


-  3  3  - 

Direction  Generation  Define  the  direction 


(4.5) 


d':  =  p(z') 


1  Phi.2')  ~  zh 
\pr.i2')  ~  2\, 


such  that  p(z')  satisfies 

(4.6)  p/,(j')=  -  uEtiIi  (>//,.  c'  +  qi,  +  L[ir1{pr1(z')  ~  z/,  )))  +  >  j  =  l,...,r 


where  lj  >  0  is  chosen  such  that  for  some  u  >  0 


(4.7) 


+Li,is)ztl  >hm2. v-' 


j  =  1,  •••,'• 


Stop  if  d'  =  0,  else  continue. 

Stepsize  Generation  z,+1  =  z'  +  A 1  cf‘ 
where 


(4.8)  f(z'  +  A V)  =  min {/(z'  +  A(f’)|ri  +  Ad’  >  0} 


4.2  Remark  The  principal  part  of  this  algorithm  consists  of  the  direction  generation  part 
(4.C),  which  can  be  performed  in  parallel  on  r  processors.  Once  this  is  done  the  stepsize 
generation  (4.8)  is  performed  and  the  new  value  z‘+1  is  shared  between  the  r  processors. 

The  following  convergence  results  were  derived  in  ( 1 3]  for  Algorithm  4.1. 

4.3  Theorem  (Convergence  of  the  Parallel  GP-SOR  Algorithm)  Let  A/  be  symmetric 
and  positive  semidefinite.  Either  the  sequence  {z‘}  generated  by  the  Parallel  GP-SOR 
Algorithm  4.1  terminates  at  a  solution  of  (3.1)  or  each  of  its  accumulation  points  solves 
(3.1). 

4.4  Corollary  (Parallel  GP-SOR  special  cases)  Condition  4.7  of  Algorithm  4.1  holds 
under  either  of  the  following  two  assumptions: 

„  2 

0  <  uj  <  mm  mm  - — - 

;=l . r  ,e/’  £„]T|A/,,|. 

l£I, 

l*  > 


(4.9) 


(4.10) 


0  <  u  <  2,  E  =  D 


and  M  is  positive  semidefinite 


-l 

Our  parallel  augmented  Lagrangian  method  consists  of  replacing  the  Serial  SOR  Al¬ 
gorithm  3.1  by  the  Parallel  GP-SOR  Algorithm  4.1  in  Algorithm  3.2  with  option  (4.10) 
for  the  choice  of  u  and  E .  We  formally  state  the  algorithm  below. 

4.5  Augmented  Lagrangian  Parallel  SOR  Algorithm 

Identical  to  Algorithm  3.2  except  that  the  Serial  SOR  Algorithm  3.1  in  Algorithm  3.2 
is  replaced  by  the  Parallel  GP-SOR  Algorithm  4.1  with  condition  (4.10)  above  in  force. 


I 


-  35  - 

5.  Computational  Results 

Our  algorithms  were  tested  on  random  linear  programs  which  were  generated  as  fol¬ 
lows.  First  the  matrix  A  of  the  linear  program  (2.1)  was  gen.rated.  Each  of  its  nonzero 
elements  was  a  random  number  generated  by  a  uniform  distribution  on  the  interval  [-100, 
100],  The  number  of  nonzero  elements  in  each  row  was  in  accordance  to  a  prescribed 
density  and  the  random  position  of  each  nonzero  element  was  determined  according  to 
a  uniform  distribution  on  the  column  indices  of  the  matrix.  Next  a  primal-dual  solution 
vector  (x,u)  was  generated  from  a  uniform  distribution  on  the  interval  [0,  10]  with  S0% 
of  the  components  being  nonzero.  Finally  the  vectors  b  and  c  were  chosen  so  that  (x,u) 
is  optimal. 

In  Figure  1  we  give  a  summary  of  our  computational  results  for  6  problems  with  the 
number  of  constraints  varying  between  25,000  and  125,000  and  the  number  of  variables 
varying  between  100,000  and  500,000.  All  tests  were  performed  on  a  MicroVax  II  with  16 
megabytes  of  memory  and  an  expanded  disk  swap  space.  We  are  not  aware  of  any  other 
linear  programming  software  that  can  handle  problems  of  the  size  that  we  have  attempted 
on  a  comparable  machine.  MINOS  [15],  a  state-of-the-art  pivotal  linear  programming 
package,  cannot  handle  any  of  the  problems  listed  in  Figure  1  because  they  are  too  big  for 
the  machine  memory  using  the  MINOS  configuration.  The  largest  problem  attempted  on 
the  MicroVax  II  with  MINOS  was  a  problem  with  5000  variables  of  20,000  constraints  and 
a  matrix  density  of  0.2%  with  about  200,000  nonzero  elements.  MINOS  was  used  with 
the  standard  partial  pricing  and  scaling  options.  After  3150  iterations  and  49  hours  54 
minutes  of  machine  time,  the  point  was  infeasible  and  the  objective  function  was  in  error 
by  59%  of  the  exact  minimum.  By  comparison  our  Algorithm  3.2  solved  the  same  problem 
in  1  hour  and  4  seconds  with  a  primal-dual  objective  function  accuracy  of  7  figures,  and 
relative  accuracy  of  not  less  than  10-9  as  defined  in  Figure  1. 

The  Parallel  SOR  Algorithm  4.5  was  implemented  on  the  Sequent  Balance  21000,  a 
multiprocessor  that  incorporates  eight  NS32032  processors  running  at  10MHz,  each  with 
a  floating  point  unit,  memory  management  and  an  8-kbyte  cache  sharing  a  global  memory 
via  a  32-bit  wide  pipelined  bus.  The  machine  has  8-Mbytes  of  physical  memory.  The 
operating  system  DYNIX,  is  a  version  of  Berkeley  4.2  bsd  unix.  The  computational  results 
are  depicted  in  Figures  2  and  3. 
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Figure  2  shows  the  total  computing  time  versus  number  of  processors  for  four  different 
densities:  d  =  1,  2,  7  and  10  percent  for  a  linear  program  with  1000  constraints  and  4000 
variables.  All  problems  were  solved  to  a  7-figure  accuracy  of  the  primal-dual  objective 
function  and  relative  accuracy  better  than  10-7  as  defined  in  Figure  1.  We  observe  that 
the  optimal  number  of  processors,  that  is  the  one  that  solves  the  problem  in  minimum 
total  time,  increases  with  density  as  expected.  This  number  is  3  for  1%  and  2%  densities. 
6  for  7%  density,  and  7  or  more  for  10%  density.  This  means  that  for  denser  problems,  a 
larger  number  of  processors  is  needed  in  order  to  arrive  at  the  shortest  solution  time.  This 
also  means  that  for  denser  problems,  the  communication  cost  does  not  become  a  dominant 
and  hence  prohibitive  factor  until  a  larger  number  of  processors  are  used. 

In  Figure  3  we  show  results  for  the  case  with  10,000  constraints  and  40,000  variables, 
with  density  of  0.35%  and  about  1,400,000  nonzero  elements.  To  our  knowledge  this  is  one 
of  the  largest  linear  programs  solved  on  this  relatively  modest  sized  multiprocessor.  One 
of  the  reasons  that  we  were  able  to  solve  larger  problems  on  the  Micro Vax  II,  is  that  the 
latter  had  twice  the  total  memory  size  of  the  Balance  21000  and  furthermore,  the  Micro  vax 
was  essentially  a  single-user  machine  dedicated  to  the  serial  SOR  algorithm.  The  optimal 
number  of  processors  for  the  low-density  case  shown  in  Figure  3  is  4.  Just  as  in  the  cases 
of  Figure  2,  the  optimal  number  of  processors  should  increase  with  problem  density. 

We  conclude  with  some  observations  on  the  speedup  efficiency  of  our  Parallel  SOR 
Algorithm  4.5.  We  define  the  speedup  efficiency  E(r)  as  the  ratio  of  the  actual  to  the 
theoretical  speedup  of  the  algorithm  using  r  processors  instead  of  1  processor,  thus 


(5.1) 


E{ r)  :  = 


T(  1) 

rT(r) 


where  T(r)  is  the  total  time  for  solving  a  given  problem  using  r  parallel  processors. 
Figure  4  shows  the  speedup  efficiencies  for  a  typical  case  of  a  linear  program  with  1000 
constraints,  4000  variables  and  2%  density.  The  reason  why  some  efficiencies  are  over 
100%  was  pointed  out  in  (13j.  The  explanation  is  that  our  Parallel  SOR  Algorithm  4.5 
changes  with  the  number  of  processors  used,  because  the  matrix  K  defined  by  (4.4)  changes 
with  the  number  of  blocks  into  which  M  is  divided.  Thus  we  are  not  -omparing  identical 
algorithms  when  we  evaluate  the  ration  T(l)/iT(r)  of  (5.1).  Nevertheless  the  express:.--',  is 
a  valid  measure  of  efficiency  in  the  sense  of  comparing  the  theoretical  reduced  time  T<\)/r 
to  the  observed  time  T(r)  for  an  algorithm  with  a  variable  K  that  depends  on  the  partition 
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of  \1.  If  the  matrix  K  is  held  fixed  for  r  =  \  and  r  >  1,  then  we  obtain  efficiencies  for 
identical  algorithms,  and  they  would  all  be  less  than  100%.  This  was  demonstrated  in 
[13],  Nevertheless  the  present  efficiencies  of  over  100%  are  indeed  very  encouraging  and 
also  give  the  additional  and  somewhat  surprising  result  that  a  serial  implementation  of 
our  r -block  Parallel  SOR  Algorithm  4.5  on  a  single  machine,  will  give  for  some  r  a  better 
computing  time  than  the  single  block  Serial  SOR  Algorithm  3.2.  For  the  specific  case  of 
Figure  4,  a  serial  implementation  of  the  r -block  parallel  SOR  with  r=2,  3,  4  and  5  will 
be  faster  than  the  single  block  serial  SOR. 
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Fig.  1.  MICROVAX  II:  Serial  SOR  Algorithm  3.2  test  results. 
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Fig.  2.  BALANCE  21000:  Total  time  for  Parallel  SOR  Algorithm  4.5  to  solve  linear  program 
versus  number  of  processors  for  various  densities  d.  (Average  of  4  randomly  generated 
cases  with  1000  constraints  and  4000  variables.) 
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Total  Time  (d  =  0  357..  m  =  10000.  n  -  -40000) 


Fig.  3.  BALANCE  21000:  Total  time  for  Parallel  SOFt  Algorithm  4.5  to  solve  linear  program 
versus  number  of  processors.  (d=density,  m=  number  of  constraints,  n=  number  of 
variables.) 
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Density  d  =  2%, 

m  -  1000,  n  =  4000 

No. 

Time  Sec. 

Speedup 

Processes 

r 

T(r) 

Efficiency 

E(r)  =  T(l)/rT(r) 

1 

3946 

— 

2 

709 

278% 

3 

638 

206% 

4 

716 

138% 

5 

711 

111% 

6 

840 

78% 

7 

850 

66% 

Fig.  4.  BALANCE  21000;  Speedup  efficiency  E(r)  for  the  Parallel  SOR  Algorithm  4.5  for 
an  LP  with  1000  constraints  and  4000  variables  and  2%  density. 
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